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Abstract 



We investigate lamellar three-phase patterns that form during the directional solidification of 
ternary eutectic alloys in thin samples. A distinctive feature of this system is that many different 
geometric arrangements of the three phases are possible, contrary to the widely studied two-phase 
patterns in binary eutectics. Here, we first analyze the case of stable lamellar coupled growth 
of a symmetric model ternary eutectic alloy, using a Jackson-Hunt type calculation in thin film 
morphology, for arbitrary configurations, and derive expressions for the front undercooling as a 
function of velocity and spacing. Next, we carry out phase-field simulations to test our analytic 
predictions and to study the instabilities of the simplest periodic lamellar arrays. For large spacings, 
we observe different oscillatory modes that are similar to those found previously for binary eutectics 
and that can be classified using the symmetry elements of the steady-state pattern. For small 
spacings, we observe a new instability that leads to a change in the sequence of the phases. Its 
onset can be well predicted by our analytic calculations. Finally, some preliminary phase-field 
simulations of three-dimensional growth structures are also presented. 
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I. INTRODUCTION 



Eutectic alloys are of major industrial importance because of their low melting points and 
their interesting mechanical properties. They are also interesting for physicists because of 
their ability to form a large variety of complex patterns, which makes eutectic solidification 
an excellent model system for the study of numerous nonlinear phenomena. 

In a binary eutectic alloy, two distinct solid phases co-exist with the liquid at the eutectic 
point characterized by the eutectic temperature Te and the eutectic concentration Ce. If 
the global sample concentration is close to the eutectic concentration, solidification generally 
results in composite patterns: alternating lamellae of the two solids, or rods of one solid 
immersed in a matrix of the other, grow simultaneously from the liquid. The fundamental 
understanding of this pattern-formation process was established by Jackson and Hunt (JH) 
They calculated approximate solutions for spatially periodic lamellae and rods that 
grow at constant velocity f , and established that the average front undercooling, that is, the 
difference between the average front temperature and the eutectic temperature, follows the 
relation 

AT = K,v\ + ^, (1) 

where A is the width of one lamella pair (or the distance between two rod centers), v is 
the velocity of the solidification front, and Ki and K2 are constants whose value depends 
on the volume fractions of the two solid phases and various materials parameters The 
two contributions in Eq. ([T]) arise from the redistribution of solute by diffusion through the 
liquid and the curvature of the solid-liquid interfaces, respectively. 
The front undercooling is minimal for a characteristic spacing 



The spacings found in experiments in massive samples are usually distributed in a narrow 



range around Xjh [2[. However, other spacings can be reached in directional solidification 
experiments by imposing a solidification velocity that varies with time. In this way, the 
stability of steady-state growth can be probed [3j]. In agreement with theoretical expections 



^ , steady-state growth is stable over a range of spacings that is limited by the occurrence of 
dynamic instabilities. For low spacings, a large-scale lamella (or rod) elimination instability 



is observed 



a. 



For high spacings, the type of instability that can be observed depends on 



the sample geometry. For thin samples, various oscillatory instabilities and a tilt instability 
can occur, depending on the alloy phase diagram and the sample concentration. Beyond 
the onset of these instabilities, stable tilted patterns as well as oscillatory limit cycles can 

n n 

be observed in both experiments and simulations p|, |6|. For massive samples, a zig-zag 
instability occurs for lamellar eutectics {7,8], whereas rods exhibit a shape instability [9 1. 

In summary, pattern formation in binary eutectics is fairly well understood. However, 
most materials of practical importance have more than two components. Therefore, eutectic 
solidification in multicomponent alloys has received increasing attention in recent years. 
A particularly interesting situation arises in alloy systems that exhibit a ternary eutectic 
point, at which four phases (three solids and the liquid) coexist. At such a quadruple point, 
three binary "eutectic valleys", that is, monovariant lines of three-phase coexistence, meet. 
The existence of three solid phases implies that there is a far greater variety of possible 
structures, even in thin samples. Indeed, for two solids a and /3, an array ... is 

the only possibility for a composite pattern in a thin sample; the only remaining degree of 
freedom is the spacing. With an additional 7 solid, an infinite number of distinct periodic 
cycles with different sequences of phases are possible. The simplest cycles are afi^afi^ . . . 
and a/3Q;7a/3a7 . . . and permutations. Clearly, cycles of arbitrary length, and even non- 
periodic configurations are possible. An interesting question is then which configurations, if 
any, will be favored. 

In preliminary works, the occurrence of lamellar structures has been reported in experi- 
ments in massive samples 10l-ll6|. The spatio-temporal evolution in ternary eutectic systems 
was observed in thin samples (quasi-2D experiments) in both metallic 17|] and organic sys- 



tems 



18| . In both cases, the simultaneous growth of three distinct solid phases from the 



liquid with a (a/3a7), (named ABAC in Ref. 18|) stacking was observed. Measurements in 
both cases revealed that y?v was approximately constant, in agreement with the JH scaling 
of Eq. (ED. 

On the theoretical side, models that extend the JH analysis from binary to ternary 
eutectics for three different growth morphologies (rods and hexagon, lamellar, and semi- 
regular brick structures) were proposed by Himemiya et al. 19| . The relation between front 
undercooling and spacing is still of the form given by Eq. ([1]), with constants Ki and K2 that 
depend on the morphology. The differences between the minimal undercoolings for different 
morphologies were found to be small. No direct comparison to experiments was given. 
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ternary eutectic growth has also been investigated by phase-field methods in 



21| . who have studied different stacking sequences formed by a = Ag2Al, (3 = {a 



Al) and 7 = AI2CU in the ternary system Al-Cu-Ag, while transients in the ternary eutectic 
solidification of a transparent In-Bi-Sn alloy were studied both by phase field modeling and 
experiments |l7 |. 

The purpose of the present paper is to carry out a more systematic investigation of 
lamellar ternary eutectic growth. The main questions we wish to address are (i) can an 
extension of the JH theory adequately describe the properties of ternary lamellar arrays 
and reveal the differences between cycles of different stacking sequences, and (ii) what are 
the instabilities that can occur in such patterns. To answer these questions, we develop a 
generalization of the JH theory to ternary eutectics which is capable of describing the front 
undercoolings of periodic lamellar arrays with arbitrary stacking sequence. Its predictions 
are systematically compared to phase-field simulations. We use a generic thermodynamically 



consistent phase-field model [22 



While this model is known to exhibit several thin- 



interface effects which limit its accuracy 24j-l28|. we show here that we can obtain a very 
satisfying agreement between theory and simulations if the solid-liquid interfacial free energy 
is evaluated numerically. In particular, the minimum-undercooling spacings are accurately 
reproduced for all stacking sequences that we have simulated. 

The model is then used to systematically investigate the instabilities of lamellar arrays, 
in particular for large spacings. We find that, as for binary eutectics, the symmetry elements 
of the steady-state array determine the possible instability modes. Whereas the calculation 
of a complete stability diagram is not feasible due to the large number of independent 
parameters, we find and characterize several new instability modes. Besides these oscillatory 
modes that are direct analogs of the ones observed in binary eutectics, we also find a new 
type of instability which occurs at small spacings: cycles in which the same phase appears 
more than once can undergo an instability during which one of these lamellae is eliminated; 
the system therefore transits to a different (simpler) cycle. Furthermore, we also find that 
the occurrence of this type of instability can be well predicted by our generalized JH theory. 

The remainder of the paper is organized as follows. In Sec. 2, we develop the generalized 
JH theory for ternary eutectics and calculate the undercooling-spacing relationships for 
several simple cycles. In Sec. 3, the phase- field model is outlined and its parameters are 
related to the ones of the theory. Sec. 4 presents the simulation results concerning both 
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steady-state growth and its instabilities. In Sec. 5, we briefly discuss questions related to 
pattern selection and present some preliminary simulations in three dimensions. Sec. 6 
concludes the paper. 

II. THEORY 

We consider a ternary alloy system consisting of components A, B and C, which can form 
three solid phases a, /3, and 7 upon solidiflcation from the liquid /. The concentrations of 
the components (in molar fractions) are denoted by c^, cb and cc and fulflU the constraint 

ca + cb + cc = 1. (3) 

This obviously implies that there are only two independent concentration flelds. 

As is customary, isothermal sections of the ternary phase diagram can be conveniently 
displayed in the Gibbs simplex. We are interested in alloy systems that exhibit a ternary 
eutectic point: four-phase coexistence between three sohds and the liquid. The isother- 
mal cross-section at the ternary eutectic temperature is displayed in Figure [H here for the 
particular example of a completely symmetric phase diagram. 

The concentration of the liquid is located in the center of the simplex (c^ = cb = cc = 
1/3), and the three solid phases are located at the corners of the eutectic tie triangle. For 
higher temperatures, no four-phase coexistence is possible, but each pair of solid phases can 
coexist with the liquid (three-phase coexistence). Each of these three-phase equilibria is a 
eutectic, and the loci of the liquid concentrations at three-phase coexistence as a function 
of temperature form three "eutectic valleys" that meet at the ternary eutectic point. On 
each of the sides of the simplex (with the temperature as additional axis), a binary eutectic 
phase diagram is found. 

The key point for the following analysis is the temperature of solid-liquid interfaces, which 
depends on the liquid concentration, the interface curvature, and the interface velocity. 
The dependence on the concentration is described by the liquidus surface, which is a two- 
dimensional surface over the Gibbs simplex. This surface can hence be characterized by two 
independent liquidus slopes at each point. For each phase v {v = a,/3,7), we choose the 
two liquidus slopes with respect to the minority components. Thus, for the a phase, the 
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Eutectic Tie-Triangle 
P +Y+ I 



Figure 1: (Color Online) Projection of the ternary phase diagram for a model symmetric ternary 
eutectic system on the Gibbs simplex. The triangle at the center is the tie-triangle at the eutectic 
temperature where four phases a, /3, 7, and / are in equilibrium. The diagram also contains the 
information on three-phase equilibria. The liquidus lines corresponding to each of these equilibria 
("eutectic valleys") are shown by dotted lines which meet at the center of the simplex, which is also 
the concentration of the liquid at which all the three solid phases and the liquid are at equilibrium. 



interface temperature is given by the generalized Gibbs-Thomson relation, 

Tj^, - Te = mS(cB - c|) + mg(cc - 4) - T^k - ^, (4) 

where cb and cc are the concentrations in the liquid adjacent to the interface, cf and 
their values at the ternary eutectic point, and m% = ~ 

^ cc =const cs=const 

the liquidus slopes taken at the ternary eutectic point. Furthermore, Fq, = ^aiTE/La is 
the Gibbs-Thomson coefficient, with 7^; the solid-liquid surface tension and La the latent 
heat of fusion per unit volume, and /i^^^ is the mobility of the a-liquid interface. For the 
typical (slow) growth velocities that can be attained in directional solidification experiments, 
the last term, which represents the kinetic undercooling of the interface, is very small. It 
will therefore be neglected in the following. The expression for the other solid phases are 
obtained by cyclic permutation of the indices. 

In the spirit of the original Jackson-Hunt analysis, for the calculation of the diffusion field 
in the liquid, the concentration differences between solid and liquid phases are assumed to 
be constant and equal to their values at the ternary eutectic point. Since we are interested 
in ternary coupled growth, which will take place at temperatures close to Te, this should be 



a good approximation. Thus, we define 

AcJ = 4 - with j = A,B,C and iy^a,j3, 7. 

In this approximation, the Stefan condition at a u-l interface, which expresses mass conser- 
vation upon solidification, reads 

d^c, = -|Ac,^ (5) 

where 9„Cj denotes the partial derivative of Cj in the direction normal to the interface, Vn 
is the normal velocity of the interface (positive for a growing solid), and D is the chemical 
diffusion coefficient, for simplicity assumed to be equal for all the components. 

We consider a general periodic lamellar array with M repeating units consisting of phases 
{1/0,1/1,1/2, ... , i/M-i) where each i/i represents the name of one solid phase (a, /3, 7) in the 
sequence, with a repeat distance (lamellar spacing) A. The width of the j-th single sohd 
phase region is {xj — Xj^i) A, with xq — and xm — 1, and the sum of all the widths 
corresponding to any given phase is its volume fraction rji,. The eutectic front is assumed to 
grow in the z direction with a constant velocity v. 




\ \ ^ ^ \ \ \ ^ 

X x,A x.,A x,A x„A x,A x.,A x,A x^A 

qAi23 01234 

Figure 2: Two examples for periodic lamellar arrays with M = 3 and M = 4 units. 



A. Concentration fields 



First, we consider the diffusion fields of the components A,B,C ahead of a growing 
eutectic front. For the calculation of the concentration fields, the front is supposed to be 
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planar, as in the sketches of Figure [2l We make the following Fourier series expansion for 

ca and cb 

oo 

cx= X^e^'^"^-""" + c^, X = A, 5. (6) 

n=— oo 

The third concentration Cc follows from the constraint of Eq. ([3]). In Eq. kn = l-nnjX 
are wave numbers and g„ can be determined from the solutions of the stationary diffusion 
equation 

vd^cx + B^'^cx = 0, 

which yields 



For all the modes n 7^ 0, we thus have g„ ~ for small Peclet number Pe = A/^ <^ 1 with 
i = 2D/vn, which will always be the case for slow growth. The mode n = describes the 
concentration boundary layer which is present at off-eutectic concentrations, and which has 
a characteristic length scale of i. 

To determine the coefficients X„ in the above Fourier series, we assume the eutectic front 
to be at the z = position. Using the Stefan condition in Eq. ([5]) and taking the derivative 
of cx with respect to the 2;-coordinate 



dzCx\z=0 = ^ -qnXnC 



n=— 00 



integration across one lamella period A of arbitrary partitioning of phases gives 



QnXJ^^X = I 5Z e-^'^-^Ac^ dx, (7) 



3=0 "-^J 

SO that the coefficients X„,n G IN in the series ansatz, Eq. ^ follow 



Xn = J2 Ac^e-^'="^(^^+^+^^)/2 sin(fc„A(x,+i + x,)/2). (8) 

j=0 



Applying symmetry arguments for the sinus and cosinus functions, we can formulate real 
combinations of these coefficients if we additionally take the negative summation indices 
into account. We obtain 

8 ^"^ 

Xn + X_n = . . ^ Ac^ cos(A;„A(xj+i + Xj)/2) sin(/c„A(a;j+i + Xj)/2), 
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i{Xn - = - — — - ^ Ac^ sin(A;„A(xj+i + Xj)/2) sm(A;„,A(xj+i - Xj)/2). 
Herewith, Eq. (|6]) reads 

M-l oo g 

cx = Cx + Xq + - — — cos(fc„A(xj+i + Xj)/2) sin(A;„A(xj+i + Xj)/2) cos(A;„x) 

j=0 n=l 

AI-l oo g 

+ ^ ^ - — — - sin(fc„A(xj+i + a;j)/2) sm{kn\{xj+i - Xj)/2) sin^k^x). (9) 
j=0 n=l '^^^^'^^ 

The general expression for the mean concentration {cx)m ahead of the m-th phase of the 
phase sequence can be calculated to yield 

^ oo A/— 1 -|^g 

= + + ^ 5Z 5Z { \2t.2/> ^^-y sin[7rn(x^+i - x^)] x 

Xm+l X^ ._ I- A fCj^ZQYi 

n — i J — u 

X sin[7rn(xj+i — Xj)] cos[7m{xm+i + Xm ~ a^j+i — (10) 

For a repetitive appearance of a phase u in the phase sequence, the mean concentration of 
component X ahead of this phase follows by taking the weighted average of all the lamellae 
of phase z/. 



T,Z=o{cx)m{Xm+l - Xm)S^^^ _ J for Z/ = I/^ 

— — jgriT witn d^^^ ~ \ , , 



B. Average front temperature 

The average front temperature is now found by taking the average of the Gibbs-Thomson 
equation along the front, separately for each phase {a,/3 and 7): 

AT, = Te-T, = -m''s{{cB)u - c|) - m^((cc). - c^) + T^in),, (12) 

for u = a, /3, 7. Here, {k)i, is the average curvature of the solid- liquid interface which can be 
evaluated by exact geometric relations to be 

/ \ _ X/m=0 {^)miXm+l ~ Xm)Sumi' 
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and 

{k) 



sin + sin 



(•^m+l ■^m)^ 

Here, Oy^y^_^ are the contact angles that are obtained by applying Young's law at the 
trijunction points. More precisely, Oy^i,^^^^^ is the angle, at the triple point (identified by the 
intersection of the two solid-liquid interfaces and the solid-solid one), between the tangent to 
the Vrn — I interface and the horizontal (the x direction). For a triple point with the phases 
I'rn.T-'m+i and liquid, the two contact angles Oy^y^j^^.O^^j^^y^ satisfy the following relations, 
obtained from Young's law, 

cos(0^,„^^^J cos(^^„_,,^„) sin(^^„^^^, + e^„^,^„)' 

Note that, in general, 9^^^^+, i- Qu^^^u^- 

A short digression is in order to motivate the closure of our system of equations. Although 
we have not given the explicit expressions, the coefficients and i?o can be simply calculated 
by using Eq. ([7]) with n = 0. However, to carry out this calculation, the width of each 
lamella has to be given. If these widths are chosen consistent with the lever rule, that is, the 
cumulated lamellar width of phase v corresponds to the nominal volume fraction of phase v 
for the given sample concentration c^, Cg , and c^, the use of Eq. ([7]) yields Xq = c\ — Cx 
{X = A,B,C). However, this result is incorrect: the concentrations of the solids are not 
equal to the equilibrium concentrations at the eutectic temperature because solidification 
takes place at a temperature below Te. Therefore, the true volume fractions depend on 
the solidification conditions. Their determination would require a self-consistent calculation 
which is exceedingly difficult. Therefore, we will take the same path as Jackson and Hunt 
in their original paper we will assume that the volume fractions of the three phases 
are fixed by the lever rule at the eutectic temperature, but we will treat the amplitudes of 
the two boundary layers, Ao and Bq, as unknowns. As in Ref. one can expect that the 
difference to the true solution is of order Pe and therefore small for slow solidification. 

With this assumption, the equations developed above can now be used in two ways. For 
isothermal solidification, the temperatures of all interfaces must be equal to the externally set 
temperature, and the three equations AT^ = AT for u = a,/3,7, can be used to determine 
the three unknowns Aq, Bq and the velocity v of the solid-liquid front. All of these quantities 
will be a function of the lamellar spacing A. In directional solidification, the growth velocity 
in steady state is fixed and equal to the speed with which the sample is pulled from a 
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hot to a cold region. The third unknown is now the total front undercooling. In the classic 
Jackson- Hunt theory for binary eutectics, the system of equations is closed by the hypothesis 
that the average undercoolings of the two phases are equal. This is only an approximation 
which is quite accurate for eutectics with comparable volume fractions of the two solids, 
but becomes increasingly inaccurate when the volume fractions are asymmetric We will 
use the same approximation for the ternary case here, and set ATq, = AT^ = ATy = AT. 
This then leads to expressions for AT as a function of the growth speed v and the lamellar 
spacing A. 



C. Examples 



1. Binary systems 



As a benchmark for both our calculations and simulations, we consider binary eutectic 
systems with components A and B and with three phases: and liquid. 




Figure 3: Sketch of a lamellar structure in a binary eutectic system with period length M = 2. Vi 
denotes a phase in the sequence (a/3). 



Setting xq = 0,xi = rja, X2 = 1, and applying Eq. (fTOj) gives 

-| CO ^ „ 

{cxU = c^ + ^o + -J2 { {^^x - A4) sin^(vrnr^.)} 

2A 

= + Xo + —V{r]^)Acx and 
?7„£ 



2A 



V{1 - r7,)Acx 



(14) 

(15) 
(16) 
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with kn = 27m/ X, Qn ~ fcn, A/£ ^ 1, Acx = Ac^ — Ac^, and the dimensionless function 

oo ^ 



n=l ^ ' 

which has the properties Virf) = V{1 — f]) = Vij] — 1). 
Furthermore, Eq. f lT2|) together with £ = IDjv leads to 

AT„ = -m^Eo - ^V{va)m%AcB + (18) 
AT^ = -m^Ao - Ac^ + Tfs {k)p, (19) 

where (k)^ = 2 sin6'a/3/(77aA) and (k)^ = 2sin6'^Q,/(r7^A). In addition, for a binary alloy 
i?o = —^0- The unknown Aq and the global front undercooling are determined using the 
assumption of equal interface undercoolings, AT^ = AT^j. The result is identical to the one 
of the Jackson-Hunt analysis. 



2. Ternary Systems 

Next, we study ternary systems with three components {A, B, C) and four phases (a, /3, 7 
and hquid). We start with the configuration (a/37a/37. . .), sketched in Figure HI 

Cb ... ^^^^ ...... 



\u = a 













XgA x^A x^A X3A 

Figure 4: Sketch of a ternary stacking order (a/37) with period length M = 3. 
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We set xq = 0, xi = r]a, X2 = rja + r]i3 = 1 — r]y and X3 = 1 and apply Eq. (fTOi) . This yields 

{cx)a = c^ + Xo + — P(r7„)A4 + Q{r]^,r]is) A4 + Q{r]^,r]^)Acl) (20) 

2A / \ 

(cx)/3 = + Xo + — (Q(r/^, r/J + V{vp)A4 + Q{r^p, r/,) A4 j (21) 



2A / \ 

(cx), = c^ + Xo + —[Q{v,,Vo.)Ac'^x + Qiv,,Vp)^4+nV,)^ciy (22) 

Here, we have used X = A, B,C and "P is the function defined in Eq. ( |T7I) . and 

CO ^ 

Qiv,^^ ,Vu,) = Y] 7 — sin(7™r7^J sin(7m?7^, ) cos[7m{r]^^ + 77^. )] (23) 

n=l ^ ^ 

P(r7i,J and Q{ri^^,'q^^) fulfill the properties P(r7j.J = -QiVf,, -VuJ and Q{i]^^,r]^^) = 



For simplicity, we now consider a completely symmetric ternary eutectic configuration: 
a completely symmetric ternary phase diagram (that is, any two phases can be exchanged 
without changing the phase diagram) and equal phase fractions rja = rjjs = r]^ = ^, which 
implies = c^- As a consequence, Xq = 0, and Eq. (!20l) simplifies to 

2A 

{cA)a -Ci = —V{Va)iAc'X - Ac^) (24) 



{cB)a -4 = [ - A4 ) (25) 



{cc)a-4 = ^^(Acg-Acjj, (26) 

for the three components. Since, in this case, all phases have the same undercooling by 
symmetry, the front undercooling is simply given by 

AT = -^V{r]a)m%AcB + (27) 

where = (sin + sin 6'q,^). The terms Ac% — Ac^ and Ac^ — Ac^^ are identical. For 
convenience, we write the preceding equation using the term we already use for the binaries 
namely Acb = Ac^ — Ac^. 

Next, we discuss again a ternary eutectic alloy with three components and four phases, 
but now for the phase cycle {a/3a'yaf3 . . .). Furthermore, we suppose that the two lamellae 
of the a phase have equal width Xrja/I. The average concentrations {cx)m are deduced from 
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x^A xX xX xA xA 

01234 

Figure 5: Schematic drawing of a ternary eutectic system with a configuration (a/^a^aP . . .) of 
periodic length M = 4. 

the general expression in Eq{TO] and read 

{cx)a = c~ + Xo + ^ [Siva, m)^^! + Q(f , + Q( t ' ^7)^4) (28) 

{cx)p = C- + X0 + — (2Q(r^^,f)A4+P(r/^)A4 + 7^(r^^,r^,)A4) (29) 

2A / \ 
(ex), = 4 + X0 + — (2Q(r/„f)A4 + 7^(r/„r/^)A4 + P(r/,)A4j, (30) 

where X = A,B, C. Furthermore, we have introduced the short notations 

00 ^ 

'^{Vu,,Vu,) = V't — sin(7mr7,,J sin(7rn?7i,J cos(7rn) (31) 

n=l ^ ' 
00 ^ 

SiVu^^Vu,) = — vTsin^(™^^./2){l + cos(7™)cos[7rn(r7^ -r^^J]}. (32) 

n=l ^ ' 

From the general formulation of the Gibbs-Thomson equation in Eq. ( |T2|) . we determine the 
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undercoolings, 



^ ^ 2 (sin + sin 9^^) 

+ ^ (2Q(r^^, f ) Ac^ + P(r7^) Ac^ + 7^(r^^, r/,) Ac^ 



+ — 



i ( Co + ^(2Q(r/^, f )AcS + V{v^)A4 + nvp,v,)^c^ 



2 sin 6 Pa 



AT, = -ml(^Ao + |^(2Q(r^„f)Ac^ + 7^(r/„r/^)A4 + P(ry,)Acl 
+ -ml (^0 + 1^ (2Q(77„ Ac^ + 7^(r/„ r/^) a4 + V{v,)Acl 



+ r,^^. (35) 

For a symmetric phase diagram (all slopes equal, = m) one can show using the as- 
sumption of equal undercooling of all phases that an expression for the global interface 
undercooling can be derived as AT = l/3(ATc + AT^ + AT,) by elimination of the con- 
stants Aq, Bq and Cq using the relation {Aq + Bq + Cq) = 0. 

D. Discussion 

A point which merits closer attention is the question which of all the possible steady-state 
configurations exhibits the lowest undercooling. Whereas the general idea that a eutectic 
system will always select the state of lowest undercooling is wrong (see Sec. |V] below), 
an information about this point constitutes nevertheless a useful starting point. Whereas 
the general solution to this problem is non-trivial, in the following we present some partial 
insights. 

Let us, for the sake of discussion, first compute the average total curvature undercooling 
ATk of an arbitrary arrangement. Consider a configuration of period M having Ma lamella 
of the a phase, lamella of the /3 phase, and lamella of the 7 phase, where the integers 
Ma, Mf), and M^ add up to M. In a system where all the solid- liquid and solid- solid surface 
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tensions are identical, the total average curvature undercooling AT^ of each phase z/ is, 

, „ 2 sin 6'M„ 

AT," = r„— ^ (36) 

,n ^ 2 sin 6* Mb 

AT,^ = r^^^^ (37) 

^ 2sin6'Mc 

= r,— — (38) 

It is remarkable that the average curvature undercooling is independent of the individual 
widths of each lamella, but depends only on the total volume fraction and the number of 
lamellae of the specific phase. Furthermore, it is quite clear from the above examples that 
the final expression for the global average interface undercooling can always be written in 
the same form as Eq. ([1]). The second term of this expression (that is, the one proportional 
to 1/A) can be computed for the case where all Gibbs-Thomson coefficients and liquidus 
slopes are equal, and reads 

K2 _ AT^ + AT^ + at; 

A ~ 3 

^2sin0 f Ma Mb MA 
= r^T^ — + — + —. 39 

3A V Va TjP J 

For the special case of a completely symmetric phase diagram and a sample at the eutectic 
composition, Eqn.( 15^ yields 

§ = r^(M. + M, + Mj, (40) 

where we have used the fact that rja = rii3 = rj^ = 1/3. Using, Ma + Mi, + Mc = M, 
K2 _ 2 sin 9 

^— = ^y^i ■ Thus, we see that the magnitude of this term per unit lamella m an 
arrangement is the same for all the possible arrangements, irrespective of the individual 
widths of the lamella and the relative positions of the lamellae in a configuration. Moreover, 
we see that for a general off-eutectic composition, choosing the number of lamellae in the 
ratio Tja : Tjp : rj^ renders the average curvature undercoolings of all the three phases equal. 
This condition is, however, relevant only for the special case of identical solid-solid and 
solid-liquid surface tensions and equal liquidus slopes of the phases. For the case when 
the solid-liquid and solid-solid surface tensions are unequal, the curvature undercooling is 
no longer independent of the arrangement of the lamella in the configuration. Hence, the 
problem of determining the minimum undercooling configuration is complex and no general 
expression regarding the number, position and widths of lamellae can be derived. 
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Another point is worth mentioning. Under the assumption that the volume fractions 
of the sohd phases are fixed by the lever rule, the width of the three lamellae in the a/37 
cycle is uniquely fixed by the alloy concentration. However, for the a/3a7 cycle, and more 
generally for any cycle with M > 3, this is not the case any more because there have to 
be at least two lamellae of the same phase in the cycle. Whereas the cumulated width of 
these lamellae is fixed by the global concentration, the width of each individual lamella is 
not. For example, in the a/3Q;7 cycle at the eutectic concentration = = = 1/3, all 
the configurations (^, 1/3, 1/3 — ^, 1/3) for < ^ < 1/3 are admissible, where the notation 
(■,■,...) is a shorthand for the list of the lamella widths Xn+i — Xn- The number ^ is an 
internal degree of freedom that can be freely chosen by the system. With our method, the 
global front undercooling can be calculated for any value of C,- For the a^aj cycle, we found 
that the configuration with equal widths of the a phases = 1/6) was the one with the 
minimum average front undercooling. This gives a strong indication that this value is stable, 
and that perturbations of ^ around this value should decay with time. Hence, the analytic 
expressions given above for the a/3a7 cycle, which are for ^ = 1/6, should be the relevant 
ones. 



III. PHASE-FIELD MODEL 



Model 



A thermodynamically consistent phase-field model is used for the present study |22l. |23|. 
The equations are derived from an entropy functional of the form 

S (e, 0) = ^ (e, c, 0) - (^ea (0, V0) + -^w (0)^ ^ dQ, (41) 

where e is the internal energy density, c = (cj)^^ is a vector of concentration variables, K 
being the number of components, and (p = {(f)a)a=i is a vector of phase-field variables, N 
being the number of phases present in the system, (f) and c fulfill the constraints 

K N 

J]q = 1 and Y,<Po. = l, (42) 

i=l a=l 

SO that these vectors always lie in — 1- and N — 1-dimensional planes, respectively. More- 
over, e is the small length scale parameter related to the interface width, s (e, c, 4>) is the bulk 
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entropy density, a {(f), V0) is the gradient entropy density and w {(p) describes the surface 
entropy potential of the system for pure capillary-force-driven problems. 
We use a multi-obstacle potential for w {(f)) of the form 

' 16 N,N,N 

— ^ (^af}4>a(Pl3 + ^ (^af}-y(pa(pp(p-f, if G X) 

W((h) = I a,/3=l a,/3,7=l (43) 

{a<l3) (q<^<7) 

OO, elsewhere 

where ^ = {0 | ^^=1 0« = 1 and (pa > 0}, (Xa/j is the surface entropy density and aap-y is 
a term added to reduce the presence of unwanted third or higher order phase at a binary 
interface (see below for details). 

The gradient entropy density a {<p, V0) can be written as 

N,N 

a (0, V0) = ^^"Z^)]^ (44) 

a,l3=l 
(a</3) 

where qa/s = (0aV0^ — ^/jV^a) is a vector normal to the af3 interface. The function ac {qap) 
describes the form of the anisotropy of the evolving phase boundary. For the present study, 
we assume isotropic interfaces, and hence ac (g^/j) = 1. Evolution equations for c and (p 



are derived from the entropy functional t 



maximization of entropy, respectively 
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irough conservation laws and phenomenological 



23| . A linearized temperature field with positive 



gradient G in the growth direction [z axis) is imposed and moved forward with a velocity w, 

T = To + G(z - vt) (45) 

where Tq is the temperature at z = at time t = 0. The evolution equations for the 
phase-field variables read 

uedtcPa = 6 (V ■ a,v<^„ (0, WcP) - a,^„ (0, V0)) - -w,^^ (0) - ^^'^'^^ - A, (46) 



where A is the Lagrange multiplier which maintains the constraint of Eq. fl42|) for </>, and 
the constant u is the relaxation time of the phase fields. Furthermore, a,V(/)Q5 (^)(i>a-> '^■xpa and 
/,(^^ indicate the derivatives of the respective entropy densities with respect to V(pa and (pa- 
The function /(c, T) in Eq. (H6ll describes the free energy density, and is related to the 
entropy density s(c, cp; T), through the relation /(c, 0; T) = e(c, 0; T) — Ts(c, 0; T), where 
e(c, 0; T) is the internal energy density. The free energy density is given by the summation 
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over all bulk free energy contributions fa{c;T) of the individual phases in the system. We 
use an ideal solution model, 

K / N 



f{c, 0; T) = 5^ Tc, \nc, + J2 ^iL V '' K {4>) , (47) 



where 

Uc; T) = f2 (tc^ In c, + c^Lf^l-Ill) (48) 
i=i ^ / 

is the free energy density of the a solid phase, and 

K 

fi{c-T) = TY,{cM{ci)) (49) 

i=l 

is the one of the liquid. The parameters L° and T° denote the latent heats and the melting 
temperatures of the i*^ component in the a phase, respectively. We choose the liquid as the 
reference state, and hence L\ = 0. 

The function ha{4>) is a weight function which we choose to be of the form ha (</>) = 
0^(3 — 20Q,). Thus, f = fa for 0^ = 1. Other interpolation functions involving other 
components of the (f) vector could also be used, but here we restrict ourselves to this simple 
choice. 

The evolution equations for the concentration fields are derived from Eq. (14T!) . 

dtc. = -V ■ (mUc, 0) + |j M,, (c, 0) V (^f^^^^^) j • (50) 

By a convenient choice of the mobilities Mij (c, (f>), self- and interdiffusion in multicompo- 
nent systems (including off-diagonal terms of the diffusion matrix) can be modelled. Here, 
however, we limit ourselves to a diagonal diffusion matrix with all individual diffusivities 
being equal, which can be achieved by choosing 

M,,(c, 0) = A(0)q (S^j - Cj) (51) 

N K 
a=l j=l 

The terms MiQ{c,<p) = Moj(c, 0) are the mobilities for the concentration current of the 
component i due to a temperature gradient. The diffusion coefficient is taken as a linear in- 
terpolation between the phases, Di[(j)) = '^a=i ^f^a, where Df is the non-dimensionalized 
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diffusion coefficient of the i component in the a phase, using the hquid diffusivity D as 
the reference, where the diffusivities of all the components in the liquid phase are assumed 
to be equal. In the simulations we assume zero diffusivity in the solid, and take the effec- 
tive diffusivity to be Di{(p) = D^(f)i. The quantity d* = aj {R/vm) is used as the reference 
length scale in the simulations, where the molar volume Vm is assumed to be independent 
of the concentration. Here, a is one of the surface entropy density parameters introduced in 
Eq. P3|) . and the surface entropies of all the phases are assumed to be equal. The reference 
time scale is chosen to be t* = d*'^ / DK The temperature scale is the eutectic temperature 
corresponding to the three phase stability regions at the three edges of the concentration 
simplex and is denoted by T* while the energy scale is given by RT*/vm- 

B. Relation to sharp-interface theory 

In order to compare our phase-field simulations to the theory outlined in Sec. 2, we 
need to relate the parameters of the phase-field model to the quantities needed as input 
for the theory. For some, this is straightforward. For example, all the parameters of the 
phase diagram (liquidus slopes, coexistence temperatures etc.) can be deduced from the free 
energy densities of Eqs. (147|) - (H9|) in the standard way. For others, the correspondence is less 
immediate. In the following, we will discuss in some detail two quantities that are crucial 
for the theory: the surface free energies and the latent heats, both needed to calculate the 
Gibbs- Thomson coefficients in Eq. (jl]). 

The surface free energy 5"q,/3 is defined as the interface excess of the thermodynamic 
potential density that is equal in two coexisting phases. For alloys, this is not the free 
energy, but the grand potential. Indeed, the equilibrium between two phases is given by 
K conditions for K components: K — 1 chemical potentials (because of the constraint of 
Eq. (E]), only K—1 chemical potentials are independent) as well as f—J2f=[^ ^i^i^ which is the 
grand potential, have to be equal in both phases. This is the mathematical expression of the 
common tangent construction for binary alloys and the common tangent plane construction 
for ternary alloys. 

The grand potential excess has several contributions. Since / = e — Ts, we need to 
consider the entropy excess. Both the gradient term in the phase fields and the potential 
?/;(</)) present in the entropy functional give a contribution inside the interface. If, along an 
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a/3 interface, all the other phase fields remain exactly equal to zero, then this contribution 
can be calculated analytically. However, this is generally not the case: in the interface, the 
phase fields (pi^, u ^ a, P can be different from zero, which corresponds to an "adsorption" of 
the other phases. Since the grand potential excess has to be calculated along the equilibrium 
profile of the fields, the presence of extra phases modifies the value of aai3- The three-phase 
terms proportional to cTq,^^ have been included in the potential function to reduce (or even 
eliminate) the additional phases. However, the total removal of these phases requires to 
choose high values of cTq^^. Such high values (>10 times the binary constant aa/s) cause 
the interface to become steeper near the regions of triple points and lines in 2D and 3D, 
respectively, which is a natural consequence of the fact that the higher order term affects 
only the points inside the phase-field simplex where three phases are present. The thinning 
of the interfaces leads to undesirable lattice pinning, which could only be circumvented by 
a finer discretization. This, however, would lead to a large increase of the computation 
times. Therefore, if computations are to remain feasible, we have to accept the presence of 
additional phases in the interfaces. 

Furthermore, there is also a contribution due to the chemical part of the free energy 



functional. This contribution, identified for the first time in Ref. [29|, arises from the fact 
that the concentrations inside the interface (which are fixed by the condition of constant 
chemical potentials) do not, in general, follow the common tangent plane, as illustrated 
schematically in Figure |6l 

Therefore, there is a contribution to the surface free energy which is given by the following 
expressions. For binary eutectic systems (A^ = 3 phases, <p = (0q, 0^, (pi); K = 2 components 
c = {caiCb))-, the vector c is one-dimensional and we define the concentration (c^) to be 
the independent field c = (c^, 1 — ca)- Then, we have 

A/,he„. (0, c; T) = / (0, c; T) - - /x^(T) [ca - c\) , (53) 

where /i^ (T) = ^ — ^ — ^ is the chemical potential of component A. For ternary eutectic 

Ca 

systems (A^ = 4 phases, = {(pa, (pi3, 07, (pi)', K = 3 components, c = (ca, cb, cc)), the vector 
c is two-dimensional and with the concentrations of A, B as the independent concentration 
fields, we get c = (c^, c^, 1 — ca — cb) and the chemical free energy excess becomes 



22 



/ 




Figure 6: (Color Online) Illustration of the existence of an excess interface energy contribution 
from the chemical free energy. Upper panel: the concentration inside the interfacial region does 
not necessarily follow the common tangent line. Here, the two convex curves are the free energy 
densities of the individual phases in contact, the straight line is the common tangent, and the 
thick non-monotonous line is the concentration along a cut through the interface. Lower panel: 
the grand chemical potential in the interface differs from the one obtained by a weighted sum of 
the bulk phase free energies, where the weighting coefficients are the interpolating functions of the 
order parameters. 



A/,,,^ (0, c; T) = /(0,c;T) - - (/.^(T)) (c^ - c^^) - (/.^ (T)) (c^ - 4) . (54) 

The entire surface excess can thus be written as the following 

5az = / (Tea {(f), V0) + -w (0) + A/^hem (0, c; T) ) dx (55) 

where x is the coordinate normal to the interface, and the integral is taken along the equi- 
librium profile (pix), c{x). This integral cannot be calculated analytically. Therefore, we 
determine the surface free energy numerically. To this end, we perform one- dimensional sim- 
ulations to determine the equilibrium profiles of concentration and phase fields, and insert 
the solution into the above formula to calculate a. For these simulations, the known bulk 
values of the concentration fields are used as boundary conditions. To accurately calculate 
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the surface excesses, it is important to include the contribution of the adsorbed phases. For 
this, the above calculations are performed by letting a small amount of these phases equili- 
brate at the interface of the major phases. Since the adsorbed phases equilibrate with very 
different concentrations compared to that of the bulk phases, the domain is chosen large 
enough such that the chemical potential change of the bulk phases during equilibration is 
kept negligibly low. 

Another important quantity which is required as an input in the theoretical expressions 
is the latent heat of fusion of the a phase. We follow the thermodynamic definition for 
the latent heat of transformation L° , 



L° = Te [s' - , (56) 
dT 



K 



and in particular = c[ln (c[) (58) 



i=l 
K 



and = J2c?^ + crin{cr), (59) 



=1 » 



where the concentrations of the phases are taken from the phase diagram at the eutectic 
temperature. 

Finally, let us give a few comments on the interface mobility yUjnt that appears in Eq. (j4]). 



In early works [30| , it was shown that an expression for this mobility in terms of the phase- 
field parameters can be easily derived in the sharp-interface limit in which the interface 
thickness tends to zero. Later on. Karma and Rappel (2j| proposed the thin-interface limit, 
in which the interface width remains finite, but much smaller than the mesoscopic diffusion 
length of the problem. This limit relaxes some of the stringent requirements of the sharp- 
interface method for the achievement of quantitative simulations. Additionally, this method 
introduces a correction term to the original expression for the interface mobility, which makes 
it possible to carry out simulations in the vanishing interface kinetics (infinite interface 
mobility) regime. 

Clearly, such modifications of the interface kinetics are also present in our model, where 
they arise both from the presence of adsorbed phases in the interface and from the structure 
of the concentration profile through the interface. Furthermore, it is well known that solute 
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trapping also occurs in phase- field models of the type used here [3l|. Since the interface 
profile can only be evaluated numerically, and since several phase-field and concentration 
variables need to be taken into account, it is not possible to evaluate quantitatively the 
contribution of these effects to the interface mobility. However, this lack of knowledge does 
not decisively impair the present study since we are mainly interested in undercooling versus 
spacing curves at a fixed interface velocity. At constant velocity, the absolute value of the 
interface undercooling contains an unknown contribution from the interface kinetics, but 
the relative comparison between steady states of different spacings remains meaningful. In 
addition, even though our simulation parameters correspond to higher growth velocities than 
typical experiments, it will be seen below that the value of the kinetic undercooling in our 
simulations is small. This indicates once more that our comparisons remain consistent. 



IV. SIMULATION RESULTS 



In this section, we compare data extracted from phase-field simulations with the theory 
developed in SecHIl for the case of coupled growth of the solid phases in directional solid- 
ification. The simulation setup is sketched in Figure [3 Periodic boundary conditions are 
used in the transverse direction, while no-flux boundary conditions are used in the growth 
direction. The box width in the transverse direction directly controls the spacing A. The 
box length in the growth direction is chosen several times larger than the diffusion length. 
The diffusivity in the solid is assumed to be zero. A non-dimensional temperature gradient, 
G is imposed in the growth direction and moved with a velocity f , such that the temperature 
field is given by Eq. (145|) . 

The outline of this section is as follows: first, we will briefly sketch how we extract the 
front undercooling from the simulation data. Then, this procedure will be validated by 
comparisons of the results to analj^ically known solutions as well as to data for binary 
alloys, for which well-established benchmark results exist. We start the presentation of our 
results on ternary eutectics by a detailed discussion of the two simplest possible cycles, a/37 
and a/3a7. We compare the data for undercooling as a function of spacing to our analytical 
predictions and determine the relevant instabilities that limit the range of stable spacings. 
Finally, we also discuss the behavior of more complicated cycles, for sequences up to length 
M = 6. 
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No Flux Boundary Conditions 




Figure 7: Simulation setup for the phase-field simulations of binary and ternary eutectic systems. 
We impose a temperature gradient G along the z direction and move it with a fixed velocity. The 
average interface position follows the isotherms at steady state in case of stable lamellar coupled 
growth. 

A. Data extraction 

At steady state, the interface velocity matches the velocity of the isotherms. The un- 
dercooling of the solid-liquid interface is extracted at this stage by the following procedure. 
First, a vertical line of grid points is scanned until the interface is located. Then, the precise 
position of the interface is determined as the position of the level line (pa = <Pi3 for an Q-P- 
interface (and in an analogous way for all the other interfaces). This is done by calculating 
the intersection of the phase-field profiles of the corresponding phases, which are extrap- 
olated to subgrid accuracy by polynomial fits. In the presence of adsorbed phases at the 
interface, the two major phases along the scan line are used for determining the interface 
point. The major phases are determined from the maximum values that a particular order 
parameter assumes along the scan line. The temperature at a calculated interface point is 
then given by Eq. ( H5l) . 

In order to test both our data extraction methods and our calculations of the surface 
tensions, we have performed the following consistency check. For an alloy with a symmetric 
phase diagram at the eutectic concentration, a lamellar front has an equilibrium position 
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when a small temperature gradient {G = 0.001) is applied to the system at zero growth 
speed. Since the concentration in the liquid is uniform for a motionless front, according to 
the Gibbs-Thomson relation the interface shapes should just be arcs of circles. This was 
indeed the case in our simulations, and the fit of the interface shapes with circles has allowed 
us to obtain the interface curvature and the contact angles with very good precision. The 
extraction of the data is illustrated in Figure [81 

120 F ' I ' ' I' - ' ^ 




Coordinate in Transverse Direction 

Figure 8: (Color Online) Procedure to extract the interface points from the simulation data with 
sub-grid resolution using higher order interpolation of the phase-field profiles. For the evaluation 
of the equilibrium properties, the solid-liquid interface points of each lamella are fitted with a circle 
which is then used to measure the radius of curvature of the particular lamella. We also calculate 
the triple point angles as the angles between the tangents to the circles at one of the points of 
intersection. 



We fit the radius and the coordinates of the circle centers. Then, the angle at the trijunc- 
tion point 6 is deduced from geometrical relations, with d = a + b and a 



Rl- Rb+ 



d — a 



e = cos"^ ( ^1 + cos"^ I ^ 



The meaning of the lengths a and h is given in Figure [HI 



2d 
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B. Validation: Binary Systems 



For comparison with the AT — A relationship known from Jackson- Hunt (JH) theory, we 
create two binary eutectic systems by choosing suitable parameters Lf and Tj" in the free 
energy density f{(l),c;T). A symmetric binary eutectic system, shown in Figure 9(a), is 
created by 



A B 
a 4.0 4.0 
/3 4.0 4.0 



\ 



( 



A B 
a 1.0 0.75980 
(3 0.75980 1.0 



\ 



A-B eutectic phase diagram 



C-D eutectic phase diagram 




0.85 




0.85 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Concentration(A) 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Concentration(C) 



(a) (b) 

Figure 9: Binary eutectic phase diagrams for a model system with stable (solid lines) and 
metastable (light dashed lines) extensions of the solidus and the liquidus lines, of (a) a symmetric 
A-B and (b) an unsymmetric C-D system. 



To create an asymmetric binary eutectic system, shown in Figure 9(b), we choose 



( 



C D 
a 5.0 5.0 
13 5.0 5.0 



C D 
a 0.96 0.80137 
(3 0.76567 1.0 



The numbers Lf , T" are chosen such that the widths of each of the (lens-shaped) two- 
phase coexistence regions remain reasonably broad, and that the approximation of using 
the values of concentration difference between the solidus and liquidus (AcJ,) at the eutectic 
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temperature for the theoretical expressions holds for a good range of undercoolings. This 
implies that the value of the Lf should not be too small. Conversely, a too high value is 
also not desirable since for large values of Lf the chemical contribution to the surface free 
energy becomes large, which leads to very steep and narrow interface profiles. 



Table I: Parameters for the sharp-interface theory, with proper calculation of the surface tension 
in the phase-field simulations for (a) a symmetric binary eutectic system with components A and 
B and (b) for an unsymmetric binary eutectic system with components C and D. 



(a) 



'5al 


1.01146 




1.01146 




1.23718 




37.70 




37.70 




4.0 


LP 


4.0 


a P 

rrv^ = rrij^ 


-0.206975 



(b) 





0.97272 


'5pi 


1.07235 


O'aP 


1.24836 


Gap 


33.903 


Spa 


41.161 




4.686 




4.711 




-0.13161 


p 


-0.22138 



We perform simulations at two different velocities V = 0.01 and V = 0.02, with a mesh 
size Ax = 1.0 and the parameter set e = 4.0, = = D^j = = 1.0, cTq,^ = cXai = 
aj3i = 1.0, aap-y = 10.0. To give an idea of the order of magnitude of the corresponding 
dimensional quantities, we remark that if we assume the melting temperatures to be around 
1700K and the other values to correspond to the Ni-Cu system used in the study of Warren 



et al. [32|, the length scale d* for the case of the binary eutectic system turns out to be 
around 0.2 nm and the time scale 0.04 ns. 

The corresponding parameters for the sharp-interface theory are given in Table [B The 



comparisons between our numerical results and the analytic theory are shown in Figs. 10(a) 



and 10(b) 



Consistent differences can be observed in the undercooling values between our data and 
the predictions from JH theory for both systems. The difference in undercoolings is smaller 
at lower velocities, which hints at the presence of interface kinetics. We find indeed that 
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(b) 



Figure 10: Comparison of AT — A relations resulting from the theoretical analysis and from the 
phase-field simulations at two different velocities for systems; (a) symmetric binary eutectic system 
(A-B) and (b) unsymmetric binary eutectic system (C-D). 



when we change the relaxation constant in the phase-field evolution equation by about 50 
%, the difference between the predicted and measured undercoolings is removed for the case 
of the considered symmetric binary phase diagram. This clearly shows that the interface 
kinetics is not negligible. It seems difficult, however, to obtain a precise numerical value for 
its magnitude in the framework of the present model. 

The spacing at minimum undercooling, however, is reproduced to a good degree of ac- 
curacy (error of 5 %), while the minimum undercooling has a maximum error of 10 %. It 
should also be noted that the JH theory only is an approximation for the true front un- 
dercooling. Results obtained both with boundary integral 6| and quantitative phase-field 



methods 



33l | have shown that, whereas the prediction for the minimum undercooling spac- 



ing is excellent, errors of 10 % for the value of the undercooling itself are typical. If the JH 
curve is drawn without taking into account the additional chemical contributions to the sur- 
face tension, a completely different result is obtained, with minimum undercooling spacings 
that are largely different from the simulated ones. We can therefore conclude that we have 
captured the principal corrections. 

In addition, we have performed equilibrium measurements of the angles at the trijunction 
point and of the radius of curvature of the lamellae as described in the preceding sub-section 
(jIV Ap for the symmetric eutectic system. The contact angles differ from the ones predicted 
by Young's equilibrium conditions only by a value of 0.2 degrees. The theoretical (from the 
Gibbs-Thomson equation) and measured undercoolings differ in the third decimal, with an 
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error of 0.1 %. 



C. Ternary Systems: Parameter set 

[a b c \ 

a 1.46964038 1.0 1.0 
1.46964038 1.0 
1.0 1.46964038 / 

C \ 
0.5 
0.5 
1.5/ 

We use a symmetric ternary phase diagram. The following matrices list the parameters 
Lf ,T" in the free energy / (</>, c; T) that were used to create a symmetric ternary eutectic 
system, shown in Figure [H We perform simulations with the parameter set e = 8.0, Ax = 

1.0, = = D^(j = 1.0,Cra7 = 0"/37 = CTjlB = CF al = (J pi = 0"^/ = 1.0, CTq/S/ = = aa^yl = 

apyi = 10.0 and compare with the theoretical expressions using the input parameters listed 
in Table HI 

Table II: Input parameters for the theoretical relations for the ternary eutectic system. 



O'al — CTpi — Cr-yl 


1.194035 


= CTay = (Jj3^ 


1.430923 


(^afi — ^ Pa — ^7cj — ^a'y 


36.81 




1.33 


mf^ = m^Q 


-0.91 


13 13 
= 


-0.91 


7 7 

mjj = 


-0.91 



13 1.0 
V7 1-0 



( A B 

a 1.5 0.5 
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D. Simple cycles: steady states and oscillatory instability 



We first perform simulations to isolate the regime of stable lamellar growth for the config- 
uration a/37. For this regime, we measure the average interface undercooling and compare 
it to our theoretical predictions. The results are shown in Figure [TTJ 
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0.05 

60 80 100 120 140 160 180 200 

Figure 11: Comparison between theoretical analysis and phase-field simulations at two different 
velocities for the arrangement (a/37) of ternary eutectic solids at V = 0.005 and V = 0.01. The 
demarcation shows the regions of stable lamehar growth and the critical spacing beyond which 
we observe amplified oscillatory behavior. There is a small region named "Damped Oscillations", 
which is a region where oscillations occur but die down slowly with time. 
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The agreement in the undercoolings is much better than for the binary eutectic systems, 
with a smaller dependence of undercoolings on the velocities. Consequently, both the spacing 
at minimum undercooling (error 4 % for V=0.005 and 6 % for V=0.01) and the minimum 
undercooling (error of 1-2 %), match very well with the theoretical relationships, as shown 
in Figure [TTl The equilibrium angles at the triple point also agree with the ones predicted 
from Young's law to within an error of 0.3 degrees, while the radius of curvature matches 
that from the Gibbs-Thomson relationship with negligible error (<0.5 %). 

It should be noted that the steady lamellae remain straight, contrary to the results of 



Ref. 



2l| . where a spontaneous tilt of the lamellae with respect to the direction of the 



temperature gradient was reported. This difference is due to the different phase diagrams: 
we are using a completely symmetric phase diagram and equal surface tensions for all solid- 
liquid interfaces, whereas 2l|] uses the thermophysical data of a real alloy. 
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Next, we are interested in the stability range of three-phase coupled growth. From gen- 
eral arguments, we expect a long-wavelength lamella elimination instability (Eckhaus-type 
instability) to occur for low spacings, as in binary eutectics [t]. Here, we will focus on oscil- 
latory instabilities that occur for large spacings. It is useful to first recall a few facts known 
about binary eutectics, where all the instability modes have been classified js, 6|. Lamellar 
arrays in binary eutectics are characterized (in the absence of crystalline anisotropy) by the 
presence of two mirror symmetry planes that run in the center of each type of lamellae, as 
sketched in Figure [T2]^a). Instabilities can break certain of these symmetries while other 
symmetry elements remain intact |3J|. In binary eutectics, the oscillatory 1-A-O mode is 
characterized by an in-phase oscillation of the thickness of all a (and /3) lamellae; both 
mirror symmetry planes remain in the oscillatory pattern. In contrast, in the 2-A-O mode, 
one type of lamellae start to oscillate laterally, whereas the mirror plane in the other type 
of lamellae survives; this leads to a spatial period doubling. Finally, in the tilted pattern 
both mirror planes are lost. 





1 
i 




1 










Figure 12: (Color Online) In a periodic arrangement of lamellae, we can identify certain lines 
of symmetry, as shown in (a) for a binary eutectic. Similarly, for the case of the two simplest 
configurations, (b) a/37 and (c) afia^ in a symmetric ternary eutectic system, such planes of 
symmetry exist. While in the case of a binary eutectic, the lines are mirror symmetry axes (shown 
by dash-dotted lines), in the special case of a symmetric ternary phase diagram, one can also 
identify quasi-mirror lines (dashed lines) where we retrieve the original configuration after a spatial 
reflection and an exchange of two phases. Only quasi-mirror lines exist in the a/37 arrangement, 
which are shown in (b), while both true- and quasi-mirror planes exist in the a/3a7 arrangement 
as shown in (c). 



It is therefore important to survey the possible symmetry elements in the ternary case. At 
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first glance, there seems to be no symmetry plane in the pattern. However, for our specific 
choice of phase diagram, new symmetry elements not present in a generic phase diagram 
exist: mirror symmetry planes combined with the exchange of two phases. Consider for 
example the /3 phase in the center of Figure [T^ b): if the system is reflected at its center, 
and then the a and 7 phases are exchanged, we recover the original pattern. At the eutectic 
concentration, there are three such symmetry planes running in the center of each lamella, 
and three additional ones running along the three solid-solid interfaces. Off the eutectic 
point, two of these planes survive if any two of the three phases have equal volume fractions. 

Guided by these considerations, we can conjecture that there are two obvious possible 
instability modes, sketched in Figure [T31 




(a) (b) 

Figure 13: (Color Online) Guided by the symmetry axes in the a/37 arrangement, one can expect 
two possible oscillatory modes at off-eutectic concentrations along the eutectic groove. The oscilla- 
tions in (a), which keep all the quasi- mirror planes intact, are expected to occur at concentrations 
towards the apex of the simplex along the eutectic groove. Another possibility, shown in (b) exists 
in which no symmetry plane remains, which is expected to occur at a concentration towards the 
binary edge of the simplex. 

In the first, called mode 1 in the following, two symmetry planes survive: the width 
of one lamella oscillates, whereas the two other phases form a "composite lamella" that 
oscillates in opposition of phase; the interface in the center of this composite lamella does 
not oscillate at all and constitutes one of the symmetry planes. In the second (mode 2), 
the lateral position of one of the lamellae oscillates with time, whereas the other two phases 
oscillate in opposition of phase to form a "composite lamella" that oscillates laterally but 
keeps an almost constant width. There is no symmetry plane left in this mode. 

The stability range of the coupled growth regime of the lamellar arrangement is indicated 
in Figured!] Steady lamellar growth is stable from below the minimum undercooling spacing 
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up to a point where an oscillatory instability occurs. In the region marked "damped oscil- 
lations" , oscillatory motion of the interfaces was noticed, but died out with time. Above a 
threshold in spacing, oscillations are amplified. We monitored the modes that emerged, and 
found indeed good examples for the two theoretically expected patterns, shown in Figure 
M 





(a) (b) (c) 

Figure 14: (Color Online) Oscillatory modes in simulations for the a/?7 configuration at the off- 
eutectic concentrations c = (0.32,0.32,0.36)) in (a), and c = (0.34,0.34,0.32) in (b) and (c). The 
spacings are A = 170 in (a) and (c) and A = 165 in (b). 



Mode 1 is favored for off-eutectic concentrations in which one of the lamellae is wider than 
the two others, such as c = (0.32, 0.32, 0.36). Indeed, in that case the (unstable) steady-state 
pattern exhibits the same symmetry planes as the oscillatory pattern. This mode can also 



appear when one lamella is sma//er than the two others, see Figure 14(c) We detect mode 2 



at the eutectic concentration, see Figure 15(b), However, a "mixed mode" can also occur, in 



which no symmetry plane survives, but the three trijunctions oscillate laterally with phase 



differences that depend on the concentration and possibly on the spacing, see Figs. 14(b) 



and 15 fc' 



Let us now turn to the a^a'-f cycle. We perform simulations for two different velocities 
V = 0.01 and V = 0.005. The comparison of the measurements with the theoretical analysis 
for steady-state growth is shown in Figure [161 For the purpose of analysis, predictions 
from the theory for both arrangements {aPj and a^a'-f) are also shown. Here again, the 
minimum undercooling spacings match those of the theory to a good degree of accuracy 
(error 5%, V=0.005). However, the undercooling is lower than the one predicted by JH- 
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Figure 15: The plot shows the trace of the triple points for the af3'j arrangement. The growth 
direction is compressed in these plots with respect to the transverse direction in order to better 
visuahze the modes. We get multiple modes at the eutectic concentration for the same spacing 
A = 159, shown in (a) and (b). In (a) we get back mode 1 while (b) matches well to our predicted 
mode 2. A mixed mode (c) is obtained at an off-eutectic concentration c = (0.34,0.34,0.32), at a 
spacing A = 165, which is a combination of oscillations in both the width and lateral spacing. 



theory, with a discrepancy of 4% for the case of V = 0.005, Figure 16(a) For V = 0.01, 



Figure 16(b) , simulations were not possible for a sufficient range of A to determine the 
minimum undercooling, because the width of the narrowest lamellae became comparable to 
the interface width e ~ 8.0 before the minimum was reached. However, the general trend of 
the data follows the predictions of the theory for both velocities. This was also the case for 
simulations carried out at an off-eutectic concentration c = (0.32, 0.34, 0.34) at a velocity of 
V = 0.005, for the same configuration af3a'y. 

Concerning the oscillatory instabilities at large spacings, it is useful to consider again the 
symmetry elements. For this cycle, there are two real symmetry planes in the steady-state 
pattern that run through the centers of the (3 and 7 lamellae. Note that these symmetries 
would exist even for unsymmetric phase diagrams and unequal surface tensions. Therefore, 
by analogy with binary eutectics, one may expect oscillatory modes that simply generalize 



the 1-A-O and 2-A-O modes of binary eutectics, see Figure [17 ( a) [ Indeed, for our simulations 

as in our 



at the eutectic concentration, we retrieve the 1-A-O type oscillation, figure 18(a) 
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Figure 16: Theoretical analysis and phase-field simulations: Comparison between the arrangements 
(a/37) and (a/3a7) for two different velocities (a) V = 0.005 and (b) V = 0.01. Plots convey 
information on the stability ranges, and the onset of oscillatory behavior of the 1-A-O type. 



hypothesis (figure 17(a)). 

This oscillatory instability can be quantitatively monitored by following the lateral posi- 
tions of the solid-solid interfaces with time. More specifically, we extract the width of the P 
phase as a function of the growth distance z. This is then fitted with a damped sinusoidal 
wave of the type Ao + Aexp{—Bz) cos((27r2;/L) + D). The damping coefficient B is obtained 
from a curve fit and plotted as a function of the spacing A. The onset of the instability is 
characterized by the change in sign of the damping coefficient. 



For the off-eutectic concentration we get two modes (figure [TSi) . While (figure 18(b)) 



corresponds well to our hypothesis to the 2-A-O type oscillation (figure 17(b)), we also 
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(a) (b) 

Figure 17: (Color Online) Predictions of oscillatory modes for the a/3a7 arrangement, reminiscent 
of the 1-A-O mode (a) and 2-A-O mode (b) in binary eutectics. 




(a) (b) (c) 

Figure 18: (Color Online) Simulations of oscillatory modes of the 0^07 configuration. The modes 
in (a) and (b) show resemblance to the 1-A-O and 2-A-O oscillatory modes of binary eutectics, 
respectively. Additionally, other modes (c) can also be observed, depending on the initial condi- 
tions. While we observe (a) at the eutectic concentration, (b) and (c) are modes at off-eutectic 
concentrations c = (0.32,0.34,0.34). The spacings are (a) A = 201, (b) A = 174 (c) A = 210. 



observe another mode as shown in figure 18(c), which combines elements of the two modes: 



both the width and the lateral position of the a lamellae oscillate. 
E. Lamella elimination instability 

For the afia'y cycle, there is also a new instability, which occurs for low spacings. We find 
that all spacings below the minimum undercooling spacing, as well as some spacings above it, 
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are unstable with respect to lamella elimination: the system evolves to the a^'j arrangement 
by eliminating one of the a lamellae, both at eutectic and off-eutectic concentrations. The 



points plotted to the left of the minimum in Figure 16(b) are actually unstable steady states 
that can be reached only when the simulation is started with strictly symmetric initial 
conditions and the correct volume fractions of the solid phases. 

This instability can actually be well understood using our theoretical expressions. As al- 
ready mentioned before, when we consider the cycle a^aj at the eutectic concentration with 
a lamella width configuration (^, 1/3, 1/3 — 1/3), the global average front undercooling at- 
tains a minimum for the symmetric pattern C, = 1/6. However, the global front undercooling 
is not the most relevant information for assessing the front stability. More interesting is the 
undercooling of an individual lamella, because this can give information about its evolution. 
More precisely, consider the undercooling of one of the a lamellae as a function of C,- If 
the undercooling increases when the lamella gets thinner, then the lamella will fall further 
behind the front and will eventually be eliminated. In contrast, if the undercooling decreases 
when the lamella gets thinner, then the lamella will grow ahead of the main front and get 
larger. A similar argument has been used by Jackson and Hunt for their explanation of the 
long- wavelength elimination instability [1]. It should be pointed out that the new instability 
found here is not a long-wavelength instability, since it can occur even when only one unit 
cell of the cycle is simulated. 

Following the above arguments, we have calculated the growth temperature of the first 
a lamella as a function of ^ using the general expressions in Eq. |33l In Figure [191 we plot 
the variation of dAT/d^ at ^ = 1/6, as a function of A. The point at which dAT/d^ 
becomes positive then indicates the transition to a stable a^a'j cycle. This criterion is in 
good agreement with our simulation results. This argument can also be generalized to more 
complicated cycles (see below). 
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Figure 19: Plot of dAT/d^, taken at ^ = 1/6 versus A for the ai/3a27 cycle, where AT is the 
undercooling of the ai lamella and ^ its width (relative to A), calculated by our analytical expres- 
sions in the volume fraction configuration. (,^, 1/3, 1/3 — ^, 1/3) at V=0.01. The cycle is predicted 
to be unstable to lamella elimination if dAT/d^, < 0. The A at which dAT/d^ changes sign is 
the critical point beyond which the a/3a7 arrangement is stable with respect to a change to the 
sequence 0/^7 through a lamella elimination. 



F. Longer cycles 

Let us now discuss a few more complicated cycles. The simple cycles we have simulated 
until now were such that during stable coupled growth the widths of all the lamellae corre- 
sponding to a particular phase were the same. This changes starting from period M — 5, 
where the configuration a^aP'-f is the only possibility (up to permutations). If we con- 
sider this cycle at the eutectic concentration and note the configuration of lamella widths 
as (^,1/3 — ^,1/3 — ^,^,1/3) and compute the average front undercooling by our theoret- 
ical expressions, we find that the minimum occurs for ^ close to 0.12. In addition, for this 
configuration, the undercooling of any asymmetric configuration (permutation of widths of 
lamellae) is higher than the one considered above. If we rewrite symbolically this configu- 
ration as tti, /32, 7, it is easy to see that this configuration has two symmetry axes of 
the same kind as discussed in the preceding subsection: mirror reflection and exchange of 
the phases a and (3. One of them runs along the interface between and a2, and the other 
one in the center of the 7 lamella. 
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Not surprisingly, our simulation results confirm the importance of this symmetry. The 
volume fractions in steady-state growth are close to those that give the minimum average 
front undercooling, see Figs. 20(b) and 20(c) 




(a) (b) (c) 

Figure 20: (Color Online) Simulations at spacings A = 135 in (a), A = 150 in (b) and A = 180 
in (c), starting from an initial configuration of af3af3'y. There is no spacing for which the a^aPj 
develops into a stable lamellar growth front. Smaller spacings switch to the a(3j arrangement while 
the larger spacings exhibit oscillatory instabilities in both the width and the lateral positions of 
the lamellae. 



Additionally, we observe oscillations in the width of the largest 7 phase and oscillations 
in the widths and the lateral position of the smaller lamellae of the a and /3 phases, while the 
interface between the larger a and (3 phase remains straight, such that the combination of 
all the a and (3 lamellae oscillates in width as one "composite lamella" . Thus, the symmetry 
elements of the underlying steady state are preserved in the oscillatory state. 

For smaller spacings, this configuration is unstable, and the sequence changes to the /3aj 



arrangement as shown in Figure 20(a) by two successive lamella eliminations. It is notewor- 
thy that we did not find any unstable sequence which switches to the a^a'y arrangement, 
which again can be understood from the presence of the symmetry. Indeed, a symmetrical 
evolution would result in a change to a configuration ai/3i7 or /32a27, but precludes the 
change to a configuration of period length M = 4. 

Going on to cycles with period M = 6, the first arrangement we consider is q;i/9i 0:2/52017, 
where we name the lamellae for eventual discussion and ease in description according 
to the symmetries. Indeed, this arrangement has two exact mirror symmetry planes in 
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the center of the 02 and the 7 phases. We find that, if we calculate the average in- 
terface undercooling curves by varying the widths of individual lamella with the con- 
straint of constant volume fraction, by choosing different ^, in the width configuration 
(^, 1/6, 1/3 — 2^, 1/6,^, 1/3), the average undercooling at the growth interface is minimal 
for the configuration (1/9, 1/6, 1/9, 1/6, 1/9, 1/3). This arrangement has the highest under- 



cooling curve among the arrangements we have considered, shown in Figure 21 (a 
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(b) 

Figure 21: (a) Synopsis of the theoretical predictions for the undercooling versus spacing of possible 
arrangements between period length M = 3 to M = 6, i.e. starting from 0/^7 to a(3"fa(3^. (b) 
Same plot but with the lamellar repeat distance A scaled with the period length M. The variation 
among the arrangements is purely a result of the variation of the solutal undercooling as can be 
infered from the discussion in Sec. IIIDI 



It also has a very narrow range of stability, and we could isolate only one spacing which 
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exhibits stable growth for A = 240, Figure 22(d) Unstable arrangements near the minimum 



undercooling spacing evolve to the a/37 arrangement, Figure 22(a), while for other unstable 



configurations we obtain the arrangements in Figure 22(b) and Figure 22(c) as the stable 



growth forms corresponding to A = 150 and A = 180 respectively. 
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Figure 22: (Color Online) Simulations starting from the arrangement al3af3aj for spacings A = 135 
in a), A = 150 in b), A = 180 in c) and A = 240 in d). 

Apart from the (trivial) period-doubled arrangement aPjaPj, another possibility for 
M = 6 is aP'-fa'-fP with a volume fraction configuration (1/6, 1/6, 1/6, 1/6, 1/6, 1/6). Simu- 
lations of this arrangement show that there exists a reasonably large range of stable lamellar 
growth, and hence we could make a comparison between simulations and the theory. We 
find similar agreement between our measurements and theory as we did previously for the 
arrangements a^j and a^aj. The plot in Figure [21] shows the theoretical predictions of all 
the arrangements we have considered until now. 



G. Discussion 

It should by now have become clear that there exists a large number of distinct steady- 
state solution branches, each of which can exhibit specific instabilities. In addition, the 
stability thresholds potentially depend on a large number of parameters: the phase diagram 
data (liquidus slopes, coexistence concentration), the surface tensions (assumed identical 
here), and the sample concentration. Therefore, the calculation of a complete stability 
diagram that would generalize the one for binary eutectics of Ref. 6|] represents a formidable 
task that is outside the scope of the present paper. Nevertheless, we can deduce from our 
simulations a few guidelines that can be useful for future investigation. 

Lamellar steady-state solutions can be grouped into three classes, which respectively 
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have (I) equal number of lamellae of all three phases (such as a/37 a/37a7/3), (II) equal 
number of lamellae for two phases (such as af3a'j), and (III) different numbers of lamellae 
for each phase. 

For equal global volume fractions of each phase (as in most of our simulations), class III 
will have the narrowest stability ranges because of the simultaneous presence of very large 
and very thin lamella in the same arrangement, which make these patterns prone to both 
oscillatory and lamella elimination instabilities. 

Any cycle in which a phase appears more than once can transit to another, simpler one 
by eliminating one lamella of this phase. This lamella instability always appears for low 
spacings below a critical value of the spacing that depends on the cycle. The possibility 
of a transition, however, depends also on the symmetries of the pattern. For instance, 
the arrangement a(3a(3a'y, if unstable, can transform into the a/37, a/3a7 or the a/3a/37 
arrangements, while for an arrangement a/3a/37, it is impossible to evolve into the a/3a7 
arrangement if the symmetry of the pattern is preserved by the dynamics. 

For large spacings, oscillatory instabilities occur and can lead to the emergence of satu- 
rated oscillatory patterns of various structures. The symmetries of the steady states seem to 
determine the structure of these oscillations, but no thorough survey of all possible nonlinear 
states was carried out. 



V. SOME REMARKS ON PATTERN SELECTION 

Up to now, we have investigated various regular periodic patterns and their instabili- 
ties. The question which, if any, of these different arrangements, is favored for given growth 
conditions, is still open. From the results presented above, we can already conclude that 
this question cannot be answered solely on the basis of the undercooling-vs-spacing curves. 
Indeed, we have shown that by appropriately choosing the initial conditions, any stable 
configuration can be reached, regardless of its undercooling. This is also consistent with 
experiments and simulations on binary eutectics [3|, |9|. To get some additional insights on 
what happens in extended systems, we conducted some simulations of isothermal solidi- 
fication where the initial condition was a random lamellar arrangement. More precisely, 
we initialize a large system with lamellae of width A = 25 and choose a random sequence 
of phases such that two neighboring lamellae are of different phases as shown in Figure 
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23(a) The global probabilities of all the phases are 1/3, which corresponds to the eutectic 
concentration, and the temperature is set to T = 0.785. 




Figure 23: (Color Online) Two snapshots of 2D dynamics in a large system. Isothermal simulations 
are started from a random configuration in (a) where the probability of occurrence of each phase is 
1/3, which is also the global concentration in the liquid. The temperature of the system is T=0.785 
and the concentration of the liquid is the eutectic concentration. A slowly changing pattern with 
a non-planar front is achieved. Some lamellae are eliminated, but no new lamellae are created. 



Under isothermal growth conditions, one would expect that, at a given undercoohng, the 
arrangement with highest local velocity would be the one that is chosen. However, in order 
for the front to adopt this pattern, a rearrangement of the phase sequence is necessary. In 
our simulations, we find that lamella elimination was possible (and indeed readily occurred). 
In contrast, there is no mechanism for the creation of new lamellae in our model, since we 
did not include fluctuations that could lead to nucleation, and the model has no spinodal 
decomposition that could lead to the spontaneous formation of new lamellae, as in Ref. 35 |. 
As a result, some of the lamellae became very large in our simulations, which led to a non- 



planar growth front, as shown in Figure 23(b) No clearcut periodic pattern emerged, such 
that our results remain inconclusive. 

We believe that lamella creation is an important mechanism required for pattern adjust- 
ment. In 2D, nucleation is the only possibility for the creation of new lamellae. In contrast, 
in 3D, new lamellae can also form by branching mechanisms without nucleation events, since 
there are far more geometrical possibilities for two-phase arrangements [36|, |37[. Therefore, 
we also conducted a few preliminary simulations in 3D. 

The cross sections of the simulated systems are 150 x 150 grid points for results in Figs. 
I2^ a) and (c), and 90 x 90 grid points for the system Figure l2^ b). The longest run took 



45 



(a) (b) (c) 

Figure 24: (Color Online) Cross-sections of patterns obtained in three-dimensional directional 
solidification. In each picture, the simulation unit cell is tiled in a 4x4 array to get a better view of 
the pattern. The pattern in (a) was started from a random configuration and evolved to a perfectly 
hexagonal pattern (at the eutectic composition for a symmetric phase diagram). At an off-eutectic 
concentration, starting with two isolated rods of a and /3 phase, the result shown in (b) is one of 
the possible structures, while with an asymmetric phase diagram at the eutectic concentration, we 
get a regular brick structure (c) from a random initial condition. 



about 7 weeks on 80 processors, for the simulation of the pattern in Figure [2^ a). This long 
simulation time is due to the fact that the pattern actually takes a long time to settle down to 
a steady state; the total solidification distance was of the order of 800 grid points. The other 
simulations required less time to reach reasonably steady states. The patterns shown in Figs. 
21] (a) and (c) start from random initial conditions (very thin rods of randomly assigned 
phases), the former with the symmetric phase diagram used previously, the latter with a 
slightly asymmetric phase diagram constructed with the changed parameters listed below, 

/ A B C \ 



a 2.0 1.0 1.0 
(5 1.0 2.0 1.0 
\7 1.0 1.0 2.0 J 



( A 



B 



C 



a 1.0 0.59534 0.63461 
(3 0.59534 1.0 0.63461 
\7 0.59534 0.59534 1.0 / 

The picture of Figure [Ml(b) corresponds to a pattern resulting of a simulation which is 
started with two isolated rods of a and /3 in a matrix of 7, with an off-eutectic concentration 
of c= (0.3,0.4,0.3). 

As shown in Figure [2ll many different steady-state patterns are possible in 3D. Not 
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surprisingly, the type of pattern seen in the simulations depends on the concentration and 
on the phase diagram. Patterns very similar to Figure [2^c) have recently been observed in 
experiments in the Al-Ag-Cu ternary system [38|. It should be stressed that our pictures 
have been created by repeating the simulation cell four times in each direction in order to 
get a clearer view of the pattern. This means that in a larger system, the patterns might 
be less regular. Furthermore, we certainly have not exhausted all possible patterns. A more 
thorough investigation of the 3D patterns and their range of stability is left as a subject for 
future work. 



VI. CONCLUSION AND OUTLOOK 



In this paper, we have generalized a Jackson-Hunt analysis for arbitrary periodic lamellar 
three-phase arrays in thin samples, and used 2D phase-field simulations to test our predic- 
tions for the minimum undercooling spacings of the various arrangements. For the model 
used here the value of the interface kinetic coefficient cannot be determined, which leads 
to some incertitude on the values of the undercooling, but this does not influence our prin- 
cipal findings. When the correct values of the surface free energy (that take into account 
additional contributions coming from the chemical part of the free energy density) are used 
for the comparisons with the theory, we find good agreement for the minimum undercooling 
spacings for all cycles investigated. Moreover, we find that, as in binary eutectics, all cycles 
exhibit oscillatory instabilities for spacings larger than some critical spacing. The type of 
oscillatory modes that are possible are determined by the set of symmetry elements of the 
underlying steady state. 

We have repeatedly made use of symmetry arguments for a classification of the oscilla- 
tory modes. In certain cases, the symmetry is exact and general, which implies that the 
corresponding modes should exist for arbitrary phase diagrams and thus be observable in 
experiments. For instance, the mirror symmetry lines in the middle of the a lamellae in 
the a/3a7 arrangement exist even for non-symmetric phase diagrams and unequal surface 
tensions, and hence the corresponding oscillatory patterns and their symmetries should be 
universal. In other cases, we have used a symmetry element which is specific to the phase 
diagram used in our simulations: a mirror refiection, followed by an exchange of two phases. 
For a real alloy, this symmetry obviously can never be exactly realized because of asym- 



47 



metries in the surface tensions, mobilities, and liquidus slopes, and therefore some of the 
oscillatory modes found here might not be observable in experiments. However, their oc- 
currence cannot be completely ruled out without a detailed survey, and we expect certain 
characteristics to be quite robust. For instance, we have repeatedly observed that two neigh- 
boring lamellae of different phases can be interpreted as a "composite lamella" that exhibits 
a behavior close to the one of a single lamella in a binary eutectic pattern. Such behavior 
could appear even in the absence of special symmetries, and thus be generic. 

Furthermore, a new type of instability (absent in binary eutectics) was found, where a 
cycle transforms into a simpler one by eliminating one lamella. We interpret this instability, 
which occurs for small spacings, through a modified version of our theoretical analysis. It is 
linked to the existence of an extra degree of freedom in the pattern if a given phase appears 
more than once in the cycle. We have not determined the full stability diagram that would 
be the equivalent of the one given in Ref. {g] for binary eutectics, because of the large number 
of independent parameters involved in the ternary problem. 

We have made a few attempts to address the question of pattern selection, with incon- 
clusive results both in 2D and 3D. In 2D, the process of pattern adjustment was hindered by 
the absence of a mechanism for lamella creation, and in 3D the system sizes that could be 
attained were too small. Based on the findings for binary eutectics, however, we believe that 
there is no pattern selection in the strong sense: for given processing conditions, the patterns 
to be found may well depend on the initial conditions and/or on the history of the system. 
This implies that the arrangement with the minimum undercooling may not necessarily be 
the one that emerges spontaneously in large-scale simulations or in experiments. 

The most interesting direction of research for the future is certainly a more complete 
survey of pattern formation in 3D and a comparison to experimental data. To this end, 
either the numerical efficiency of our existing code has to be improved, or a more efficient 
model that generalizes the model of Ref. {33 1 to ternary alloys has to be developed. 
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